Discreteness and its effect on the water-wave turbulence. 
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Abstract 



in 

> 

We perform numerical simulations of the dynamical equations for free water surface in finite 
basin in presence of gravity. Wave Turbulence (WT) is a theory derived for describing statistics of 



weakly nonlinear waves in the infinite basin limit. Its formal applicability condition on the minimal 
size of the computational basin is impossible to satisfy in present numerical simulations, and the 
number of wave resonances is significantly depleted due to the wavenumber discreteness. The 



goal of this paper will be to examine which WT predictions survive in such discrete systems with 
depleted resonances and which properties arise specifically due to the discreteness effects. As in [1- 
3], our results for the wave spectrum agree with the Zakharov-Filonenko spectrum predicted within 
WT. We also go beyond finding the spectra and compute probability density function (PDF) of the 
wave amplitudes and observe an anomalously large, with respect to Gaussian, probability of strong 
waves which is consistent with recent theory [4, 5] . Using a simple model for quasi-resonances we 
predict an effect arising purely due to discreteness: existence of a threshold wave intensity above 
which turbulent cascade develops and proceeds to arbitrarily small scales. Numerically, we observe 



that the energy cascade is very "bursty" in time and is somewhat similar to sporadic sandpile 
avalanches. We explain this as a cycle: a cascade arrest due to discreteness leads to accumulation 
of energy near the forcing scale which, in turn, leads to widening of the nonlinear resonance and, 
therefore, triggering of the cascade draining the turbulence levels and returning the system to the 
beginning of the cycle. 



1 Introduction 

Normal state of the sea surface is chaotic with a lot of waves at different scales propagating in random 
directions. Such a state is referred to as Wave Turbulence (WT). Theory of WT was developed by 
finding a statistical closure based on the small nonlinearity and on the Wick splitting of the Fourier 
moments, the later procedure is often interpreted as closeness of statistics to Gaussian or/and to phase 
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randomness (the two are not the same, see [4-6]). This closure yields a wave-kinetic equation (WKE) 
for the waveaction spectrum. Such WKE for the surface waves was first derived by Hasselmann [7] . A 
significant achievement in WT theory was to realize that the most relevant states in WT are energy 
cascades through scales similar to the Kolmogorov cascades in Navier-Stokes turbulence, rather than 
thermodynamic equilibria as in the statistical theory of gases. This understanding came when Zakharov 
and Filonenko found an exact power-law solution to WKE which is similar to the famous Kolmogorov 
spectrum [8]. 

Numerical simulation of the moving water surface is a challenging problem due to a tremendous 
amount of computing power required for computing weakly nonlinear dispersive waves. This arises 
due to presence of widely separated spatial and time scales. As will be explained below, the weaker 
we take nonlinearity, the larger we should take the computational box in order to overcome the k- 
space discreteness and ignite wave resonances leading to the energy cascade through scales. In order 
to maximize the inertial range, one tries to force at the lowest wavenumbers possible, but without 
forgetting that the forcing should be strong enough for the resonance broadening to overcome the 
discreteness effect. But nonlinearity tends to grow along the energy cascade toward high fc's [9, 10] and, 
therefore, the forcing at low wavenumbers should not be too strong for the nonlinearity to remain weak 
throughout the inertial range. A simple estimate [11] says that for the resonant interactions to be fully 
efficient, one must have a N x N computational box with the number of modes N related to the mean 
surface angle a as 

N > 1/a 4 . 

Thus, for small enough nonlinearity a ~ 0.1 one must have at least 10000 x 10000 resolution which is 
far beyond present computational capacity. Thus, in none of the existing numerical experiments, nor 
in near future, the formal applicability conditions of WT can be realized. On the other hand, a small 
fraction of the resonances can be activated at levels of nonlinearity which is much less than in the above 
estimate and this reduced set of resonant modes can, in principle, be sufficient to carry the turbulent 
cascade through scales. In this paper we estimate the minimal resonance broadening which is sufficient 
to generate the cascade. 

It is necessary to examine which of the WT predictions survive beyond the formal applicability 
conditions when the cascade is carried by a depleted set of resonant modes, and which specific features 
arise due to such resonance depletion. 

As in other recent numerical experiments [1-3], here we observe formation of a spectrum consistent 
with the ZF spectrum corresponding to the direct energy cascade. For this, the resonance broadening 
at the forcing scale should be maintained at about an order of magnitude larger than the minimal level 
necessary for triggering the cascade. Further, in agreement with more recent WT predictions about the 
higher-order statistics [4, 5] , we observe an anomalously high (with respect to Gaussian) probability of 
the large-amplitude waves. Also in agreement with recent WT findings [5], we observe a buildup of 
strong correlations of the wave phases <pk whereas the factors e l ^ k remain de-correlated. 

There are also distinct features arising due to discreteness. We analyze them by exploiting the 
two-peak structure of the time-Fourier transform at each k: a dominant peak at (very close to) linear 
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frequency Uk, and a weaker one with a frequency approximately equal to lujk/2- The second peak is a 
contribution of the A;/2-mode in the nonlinear term of the canonical transformation relating the normal 
variable and the observables (e.g. surface elevation). In fact, the nature of the second frequency peak is 
quite well understood in literature and it has even been used for remote sensing of vertical sheer by VHF 
frequency radars [13]. In absence of nonlinearity one would observe only the first but not the second peak 
and, therefore, one can quantify the nonlinearity level as the ratio of the amplitudes of these two peaks. 
When the second peak becomes stronger than the first one, the wave phase experiences a rapid and 
persistent monotonic change. Detecting such phase "runs" gives an interesting picture of the nonlinear 
activity in the 2D /c-space. In particular, we notice a "bursty" nature of the energy cascade resembling 
sandpile avalanches. Possible explanation of such behavior is the following. When nonlinearity is 
weak, there is no wave-wave resonances, consequently there is no effective energy transfer, and system 
behaves like "frozen turbulence" (term introduced in [14] for the capillary wave turbulence). Energy 
generated at the forcing scale will accumulate near this scale and the nonlinearity will grow. When 
resonance broadening gets wide enough, so that the resonances are not inhibited by discreteness, the 
nonlinear wave-wave energy transfer starts, which diminishes nonlinearity and subsequently "arrests" 
resonances. Thus the system oscillates between having almost linear oscillations with stagnated energy 
and occasional avalanche- like discharges. 



2 Equations for the free surface 

Let us consider motion of a water volume of infinite depth embedded in gravity and bounded by a 
surface separating it from air at height z = ?y(x, t) where x = (x,y) is the horizontal coordinate. Let 
the velocity field be irrotational, u = V$, so that the incompressibility condition becomes 

A$ = 0, for z < ri(x,t). (1) 

Rate of change of the surface elevation must be equal to the vertical velocity of the fluid particle on 
this surface, which gives 

D t r] = d z ®, for z = ri(x,t), (2) 

where D t = d t + u • Vj_ is the material time derivative. The second condition at the surface arises from 
the Bernoulli equation in which pressure is taken equal to its atmospheric value. This condition gives 

d t $ + -|V$| 2 = -gr), for z = rj(x,t), (3) 

where g is the free-fall acceleration. 

Although equations (2) and (3) involve only two-dimensional coordinate x, the system remains 
three-dimensional due to the 3D equation (1). One can transform these equations to a truly 2D form 
by assuming that the surface deviates from its rest plane only by small angles and by truncating the 
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nonlinearity at the cubic order with respect to the small deviations. This procedure yields the following 
dynamical equations (see e.g. [15, 17]): 



r[*] 

- £ (r[r [*]iri + v_ L .[(v_ L ¥)iri) 



(4) 



+e 2 (y [r [r m V ] v] + [(A ± *) ?? 2 ] + (r [*] ^ 2 ) 



*t = -5-77 
1 



^d v ^i 2 -( r [*D 2 ) 

e 2 rm (r [r [*] v ] + (A ± *)ri), 



(5) 



where 



* = $ 



Z=7?(x,t)) 



(6) 



and T is the Gilbert transform which in the Fourier space corresponds to multiplication by k = |k|, i.e. 



In equations (4) and (5), we rescaled variables r\ and \& to make them order one, so that the nonlinearity 
smallness is now in a formal parameter e C 1. 

Truncated equations (4) and (5) will be used for our numerical simulations. They have a convenient 
form for the pseudo-spectral method which computes evolution of the Fourier modes but switches back 
to the coordinate space for computing the nonlinear terms. However, for theoretical analysis these 
equations have to be diagonalised in the /c-space and a near-identity canonical transformation must be 
applied to remove the nonlinear terms of order e since the gravity wave dispersion uj = y/gk does not 
allow three-wave resonances. The resulting equation is also truncated at e 2 order and it is called the 
Zakharov equation [16-19], 



where u l ® u = u>i + u a — uj^ — u u , uj\ = \fgk~i is the frequency of mode 1 ( 1 G Z 2 ), k/ = 2nl/L is the 
wavenumber, L is the box size and ki = \ki\. Here, a\ is the wave action variable in the interaction 
representation, a\ = e lu3lt bi where b\ is a normal variable, 




Here, we have the following convention for defining the Fourier transform 




(7) 




(8) 




(9) 
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Here, 0(e) terms appear because of the near-identity canonical transformation needed to remove the 
quadratic terms from the evolution equations [16-21]. Expression for the interaction coefficient Wj% is 
lengthy and can be found in [21]. 

Zakharov equation is of fundamental importance for theory and it is also sometimes used for nu- 
merics. However, in our work we choose to compute equations (4) and (5) because this allows us to 
use the standard trick of pseudo-spectral methods via computing the nonlinear term in the real thereby 
accelerating the code. 

3 Statistical Quantities in Wave Turbulence 

Let us consider a wavefield in a periodic square basin of side L and let the Fourier representation of 
this field be a/(t) where index \eZ 2 marks the mode with wavenumber k/ = 27rl/L on the grid in 
the 2-dimensional Fourier space. Discrete /c-space is important for formulating the statistical problem. 
For simplicity let us assume that there is a cut-off wavenumber k max so that thee is no modes with 
wavenumber components greater than k max , which is always the case in numerical simulation. In this 
case, the total number of modes is N = {k max j-nV) 2 and index 1 will only take values in a finite box, 
1 G Bn C Z 2 which is centered at and all sides of which are equal to N 1 ^ 2 . To consider homogeneous 
turbulence, the large box (i.e. continuous k) limit, N — > oo, will have to be taken later. 

Let us write the complex a/ as ai = A^i where A[ is a real positive amplitude and ipi is a phase 
factor which takes values on 5 1 , a unit circle centered at zero in the complex plane. The most general 
statistical object in WT [5] is the iV-mode joint PDF P^' defined as the probability for the wave 
intensities A 2 to be in the range (si,si + dsi) and for the phase factors ipi to be on the unit-circle 
segment between £j and + d£i for all 1 G Bn- 

The fundamental statistical property of the wavefield in WT is that all the amplitudes Ai and phase 
factors tpi are independent statistical variables and that all ^'s are uniformly distributed on S 1 . This 
kind of statistics was introduced in [4-6] and called "Random Phase and Amplitude" (RPA) field. In 
terms of the PDF, we say that the field a is of RPA type if it can be product-factorized, 

where p£ a \si) is the one-mode PDF for variable Af. 

Note that in this formulation the distributions of A\ remain unspecified and, therefore, the ampli- 
tudes do not have to be deterministic (as in earlier works using RPA) nor do they have to correspond 
to Gaussianity, 

P J (a) ( a ,) = -exp(- a ,/ni) (11) 
n i 

where rii = (si) is the waveaction spectrum. 

Importantly, RPA formulation involves independent phase factors ip = and not phases <fi. Firstly, 
the phases would not be convenient because the mean value of the phases is evolving with the rate equal 
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to the nonlinear frequency correction [5]. Thus one could not say that they are "distributed uniformly 
from —Ti to 7r" . Moreover the mean fluctuation of the phase distribution is also growing and they quickly 
spread beyond their initial 27r-wide interval [5]. But perhaps even more important, it was shown in [5] 
that 0's build mutual correlations on the nonlinear timescale whereas ^'s remain independent. In the 
present paper we are going to check this theoretical prediction numerically by directly measuring the 
properties of 0's and ^'s. 

In [4, 6] RPA was assumed to hold over the nonlinear time. In [5] this assumption was examined 
a posteriori, i.e. based on the evolution equation for the multi-point PDF. Note that only the phase 
randomness is necessary for deriving this equation, whereas both the phase and the amplitude random- 
ness are required for the WT closure for the one-point PDF or the kinetic equation for the spectrum. 
This fact allows to prove that, if valid initially, the RPA properties survive in the leading order in small 
nonlinearity and in the large-box limit [5]. Such an approximate leading-order RPA is sufficient for the 
WT closure. 



4 Theoretical WT predictions 

When the wave amplitudes are small, the nonlinearity is weak and the wave periods, determined by the 
linear dynamics, are much smaller than the characteristic time at which different wave modes exchange 
energy. In the other words, weak nonlinearity results in a timescale separation and this fact is exploited 
in WT to describe the slowly changing wave statistics by averaging over the fast linear oscillations. 



4.1 Evolution equations for the PDF's, moments and spectrum. 

In [5] the following equation for for the iV-mode PDF was obtained for the four-wave systems, 



' 5 ' 
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5s n 



(12) 
(13) 



Here N — > oo limit has already been taken and means the variational derivative. Using this equation, 
one can prove that RPA property holds over the nonlinear time, i.e. the iV-mode PDF remains of the 
product factorized form with accuracy sufficient for the WT closures to work [5]. Using RPA, we get 
for the one-point marginals [5], 

dP a dF . 

-w + arr°- (14) 



with F is a probability flux in the s-space, 

F = Sj^Pa 



Vj 



(15) 
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where 



Vj 



Aire 



[ \w jl |V 

/ I nm I w nr 



r m ^(^n m ) n i n mn n dkidk m dk n , 



Ij = 47T6 2 / \Wi? m Kj{ 
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jl 



ni{n m + n n ) - n m n n 
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Here we introduced the wave-action spectrum, 



n 3 = M 2 >. 



From (14) we get the following equation for the moments M- = (A* p ): 



Mf ] = - P1] M ( f ) + P \M ( f- l) . 
which, for p — 1 gives the standard wave kinetic equation (WKE), 



(16) 

(17) 

(18) 

(19) 
(20) 



4.2 Preservation of the RPA property. 

Validity of the WT theory relies on persistence of the RPA property of the wavefield over the nonlinear 
evolution time. Such persistence was demonstrated in [5] based on the evolution equation for the multi- 
mode PDF, where the product factorization of PDF was shown to hold with an accuracy sufficient for 
the WT closure. It was also emphasized in [5] that RPA must use independent phase factors ipk rather 
than the phases independence of which does not survive over the nonlinear time. The theoretical 
prediction of persistent independence of Ak's and ?/Vs and about the growth of correlations of 0fc's will 
be checked in this paper numerically. 

Further, the WT approach predicts that the mean value of the phase grows with a rate given by 
the nonlinear frequency correction and that the r.m.s. fluctuation of the phase also grows in time 
[5]. In this paper we will see that in reality the time evolution of the phase is more complicated than 
this WT prediction: the phase exhibits quasi-periodic fluctuations intermittent by rare "phase runs" , - 
monotonic changes over several linear periods by large values which can significantly exceed 2n. 



4.3 Steady state solutions. 

Steady power-law solutions of WKE which correspond to a direct cascade of energy and an inverse 
waveaction cascade are, 

n(Jfe) = CiP 1/3 £T 4 (21) 
n(Jfe) = C 2 Q 1/3 Ar 23 / 6 , (22) 
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where P and Q are the energy and the waveaction fluxes respectively and C\ and C2 are constants, and 
k — |k|. The first of these solutions is the famous ZF spectrum [8] and it has a great relevance to the 
small-scale part of the sea surface turbulence. It has been confirmed in a number of recent numerical 
works [1-3], but we will also confirm it in our simulation. 

Now, let us consider the steady state solutions for the one-mode PDF. Note that in the steady state 
7/77 = n which follows from WKE (20). Then, the general steady state solution to (14) is 

P = const exp (—s/n) — (F / rj) Ei(s / n) exp (— s/n), (23) 

where Ei(x) is the integral exponential function. At the tail s ^> n& we have 

F 

P^-— (24) 

57 

if F ^ 0. The 1/s tail decays much slower than the exponential (Rayleigh) part and, therefore, it 
describes strong intermittency. On the other hand, 1/s tail cannot be infinitely long because otherwise 
the PDF would not be normalizable. As it was argued in [5], the 1/s tail should with a cutoff because 
the WT description breaks down at large amplitudes s. This cutoff can be viewed as a wavebreaking 
process which does not allow wave amplitudes to exceed their critical value, P(s) = for s > s n [. 

Relation between intermittency and a finite flux in the amplitude space was observed numerically 
also for the Majda-Mc-Laughlin-Tabak model by Rumpf and Biven [22]. 



5 Resonant interaction in discrete /c-space 

Importance of the /c-space discreteness for weakly nonlinear waves in finite basins were realized by 
Kartashova [23] and it was later discussed in a number of papers, [14, 28-30]. Nonlinear wave interactions 
crucially depend on the sort of resonances the dispersion relation allows. The dispersion relation of the 
surface gravity waves is concave and hence it forbids three-wave interactions, so that the dominant 
process is four-wave. Resonant manifold is defined by resonant conditions k + ki = k2 + k3, + = 
^k 2 + ^k 3 or, substituting = y/gk, 

k + ki = k 2 + k 3 

Vk+^/h = a/A^+a/A^ (25) 

The problem of finding exact resonances in discrete /c-space was first formulated by Kartashova [23] 
who gave a detailed classification of for some types of waves, e.g. Rossby waves. For the deep-water 
gravity waves, Kartashova only considered symmetric solutions, i.e. of the form k = k 2 , £7 = k 3 (or 
k — &3, ki — £2). Interestingly, there appear to be also asymmetric solutions, 

Solutions for collinear quartets and the tridents are easy to find by rewriting the resonant conditions 
(25) as polynomial equations (appropeiately re-arranging and taking squares of the equations) and 
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using rational parametrisations of their solutions. This way we get the following family of the collinear 
quartets [24], 

k=(a,0), k 1 = (6,0), k 2 = (c,0), ks = (d,0) 

with 

a = m 2 (m + n) 2 , b = n 2 (m + n) 2 , c = —m 2 n 2 , d = (m + n) A + m 2 n 2 , 
where m and n are natural numbers, and "tridents" , 

k=(a, 0), ki = (-6, 0), k 2 = (c, d), k 3 = (c, -<Z) 

with 

a = ( s 2 + t 2 + st) 2 , b= (s 2 + t 2 - st) 2 , c = 2st(s 2 + t 2 ), d = s 4 -t 4 , 

where s and t are integers. Both of these classes can be easily extended by re-scaling, i.e. multiplying all 
four vectors by an integer. Even more solutions can be obtained via rotation by an angle with rational- 
valued cosine (this gives a new rational solution) and further re-scaling (to obtain integer solution out 
of the rational one). 

There are other, rather rare, exact nontirvial resonances, for example one provided to us by Kar- 
tashova [25]: 

k = (495, 90), ki = (64, 128), k 2 = (359, 118), k 3 = (200, 100). 

The complete investigation of the exact resonance types and their respective roles is a facinating subject 
of future work). 

Note that the interactions coefficient vanishes on the collinear quartets [27] and, therefore, these 
quartets do not contribute into the turbulence evolution. Secondly, there appears to be much more 
symmetric quartets than tridents, so the later are relatively unimportant for the nonlinear dynamics 
too. 

Because of nonlinearity, the wave resonances have a finite width. Even though this width is small in 
weakly nonlinear systems, it may be sufficient for activating new mode interactions on the fc-space grid 
and thereby trigger the turbulent cascade through scales. Such quasi-resonances can be roughly modeled 
through k x + k 2 = k 3 + k 4 , |o; kl + c<j k2 — cj k3 — cj k4 | < 5, where 5 describes the resonance broadening. 
Figure 1 shows quasi-resonant generations of modes on space [— 64, 64] 2 , where initially (generation 1) 
only modes in the ring 6 < k < 9 were present (as in our numerical experiment). With broadening 
of the resonant manifold smaller than the critical 8 cr u ~ 1.4 * 10 -5 , a finite number of modes outside 
the initial region get excited due to exact resonances (generation 2) but there are no quasi-resonances 
to carry energy to outer regions in further generations. If broadening is larger than critical, energy 
cascades infinitely. Contrary to the the case of capillary waves where quasi-resonant cascades die out 
if broadening is not large enough [28] , in the case presented here quasi-resonant cascades either do not 
happen at all (if S < 6 ait) or they spread through the wavenumber space infinitely (if 5 > 5 C rit), as 
happens on Figure (1). 
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Figure 1: Quasi-resonant generations of modes on Fourier space [— 64, 64] 2 , where initially (generation 
1) only modes in the ring 6 < k < 9 were present. Each next generation consists of the union of the 
modes of the current generation and the new modes which satisfy the quasi-resonance condition with 
the modes of the current generation. Here the value of the broadening is slightly above S crit 
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6 Numerical simulation 



Numerical simulations presented in this work were performed on a single-processor workstation (2.5GHz, 
1Gb RAM). We performed a direct numerical simulation, integrating the dynamical equations of mo- 
tion (4) and (5) using pseudo-spectral method with resolution of 256 x 256 wavenumbers. Numerical 
integrator used for advancing in time was RK7(8) presented in [31]. Time step was -^p- where T min is 
the period of the shortest wave on the axis. Approximate processor time for this work was 4.5 weeks. 

In our numerical experiment, we force the system in the /c-space ring k* < k < k* with k* = 6 
and k* = 9. This ring is located at the low wavenumber part of the /c-space in order to generate 
energy cascade toward large fc's, but we deliberately avoid forcing even longer waves (k < 6) because 
our experience shows that this would lead to undesirable strong anisotropic effects. In the ring, we 
fix the shape to coincide with the ZF spectrum, < |ak(^)| 2 >~ k~ 4 , and hence we set |?7k| ~ k~ 7 / 4 , 
l^kl ~ /c~ 9 / 4 . These fixed amplitudes were then multiplied by random phase factors. Thus, surface r?k 
and velocity potential \I/k were set to 2n 3 e ldk * k ak , where #k were uniformly distributed in [0, 2ir] and 



where x = —7/4 for the surface and x = —9/4 for velocity potential. Damping was applied in both 
small and large wavenumber regions. At small wavenumbers inside the forcing ring, we applied an 
adaptive damping to prevent formation of undesirable "condensate" which could spoil isotropy and 
locality of scale interactions. At large wavenumbers, damping is needed to absorb the energy cascade 
and, therefore, to avoid "bottleneck" spectrum accumulation near the cutoff wavenumber. In our 
simulations, we implemented the damping as a low-pass filter 7k applied to the /c-space variables at 
each time step. The damping function had the form 



Nonlinearity parameter was set to e = 2 • 10 2 , which is a sufficient value to produce a resonance 
broadening for supporting energy cascade. 

7 Results 
7.1 Spectrum 

Measuring the spectrum has by far dominated WT studies because this is the most basic and robust 
theoretical object and because this quantity is easier to observe experimentally that more subtle statis- 
tical quantities. For the surface gravity waves, the WT prediction about the energy cascade spectrum 
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Figure 2: Waveaction spectrum of waves. Dashed and solid lines show spectra after 300 and 2000 
periods of mode k = 10. The straight line corresponds to the ZF spectrum (—4 slope). 



have been confirmed in several recent numerical studies [1-3]. Here, we also start by presenting the 
spectrum. Figure 2 shows the spectrum at t — 300T P at t — 2000T P , where T p is period of the slowest 
mode (k = 10). The obtained spectrum is in agreement with the k~ A shape predicted by WT theory in 
the inertial range, and this serves as a validation of our code. Our subsequent statistical measurements 
will be made at the time where this steady state spectrum has already got established. 

7.2 Wave-amplitude probability density function and its moments 

Now we consider the PDF of amplitudes for which predictions where made recently within the WT 
approach. To measure the PDF of amplitudes |ak(t)| 2 , we set two radial regions in k-space £45 = 
[13, 17], and k^ = [33, 37]. These regions were inside the inertial range and had well mixed phases and 
amplitudes since the experiment was done after performing 2000 rotations of the peak mode. We looked 
at the time-span of approximately 855 rotations of modes k = 15 or 1230 rotations of modes k = 35 
and collected amplitudes of all modes from these three regions. The number of amplitudes collected 
was over 1.1 million in region £45 and over 2.1 million in region £35. 

Figure 3 shows a log-plot of the PDF the amplitude A\ in £45 and an exponential fit of its low- 
amplitude part. One can see intermittency, i.e., an anomalously large probability of strong waves. 
We can also see that this discrepancy from Gaussianity happens in the tail, i.e. well below the mean 
amplitude value s — n^. While the PDF tail in not long enough for drawing decisive conclusions about 
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Figure 3: Probability density function for the amplitude |afc| 2 with k G [13, 17]. The linear fit is shown 
based on the slope of the low-amplitude part (the Gaussian core). The vertical straight line marks the 
mean value (spectrum) rik = (|afe| 2 )- 
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Figure 4: Probability density function for the amplitude \dk\ 2 with k G [33,37]. Same notations as in 
Figure 3. 



realization of the theoretically predicted 1/s scaling, it certainly gives a conclusive evidence that the 
probabilities of large amplitudes are orders of magnitude higher than in Gaussian turbulence. Figure 4 
shows the PDF of A\ in k 35 . We can see some non-Gaussianity in k 35 as well, although much less than 
in £45. Similar conclusion that the gravity wave turbulence is more intermittent at low rather than high 
wavenumbers was reached on the basis of numerical simulations in [3]. 

Deviations from Gaussianity can be also seen in figure 5 which shows ratio of the moments M^> = 
(\a\ 2p ) to their values in Gaussian turbulence, np\. Again, we can see that such deviations are greater 
at the small-/c part of the inertial range. 

7.3 Frequency properties. 

We examine the frequency properties of waves by performing the time-Fourier transform at each fixed 
wavenumber. A typical plot, for k = (17,0), is shown in Figures 6. Our first observation is that we 
always see two peaks - the bigger one at the linear frequency and a smaller peak at a shifted frequency. 
We interpret the second peak as a nonlinear effect since there is no frequency shift in the linear system. 

Also, it appears that the measured ratio of squares of the peak frequencies is approximately equal 
to 2 for all wavenumbers (within 10% accuracy). This can be explained by the nonlinear term in the 
canonical transformation (9), e.g. 0(e)-term which is quadratic with respect to the wave amplitude. In 
particular, the mode k/2 makes contribution to this term which oscillates at frequency 2u;(k/2) = \fTk~ 



14 



0.8 



0.7 
0.6 
£ 0.5 

Q. 

1, 0.4 

o 




*1 2 3 4 5 6 7 8 

P 



Figure 5: Ratio of the moments = (\a\ 2p ) to their values in Gaussian turbulence, np\, for k e 
[13, 17] (dashed line), k G [23,27] (thin solid line) and k £ [33,37] (thick solid line). 
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Figure 6: Frequency distribution of waveaction at a fixed wavenumber k = (17, 0). 
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Figure 7: Graph of A\{t) at wavenumber k = (25,0). 



which appears to coincide with the second peak's frequency. Thus we see that contribution of k/2 
dominates in the nonlinear term of the canonical transformation. 

The two-frequency character at each wavenumber has an interesting relation to the amplitude and 
phase dynamics as will be seen in the next section. 

7.4 Amplitude and phase evolution. 

Figures 7 and 8 show time evolution of the amplitude Ak and the (interaction representation) phase 4>k 
respectively for k = (17,0). The phase evolution seen in figure 8 consists of the time intervals when 
it oscillates quasi-periodically (with amplitude less than 2n) which are inter-leaved by sudden "phase 
runs", - fast monotonic phase changes by values which can significantly exceed 2n. By juxtaposing 
figures 7 and 8, one can see that when the phase runs happen then the amplitude is close to zero. 
We zoom in at the amplitude graph in a time interval characterized by low amplitudes (and therefore 
phase runs) in figure 9. One can see large-amplitude quasi-periodic oscillations on A k , - it changes in 
value several-fold over a time comparable to the linear wave period. This indicates that such phase run 
intervals mark the places where the WT assumption of weak nonlinearity breaks down. 

This behavior, together with the two-peak character of the time-Fourier spectrum, suggest that at 
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Figure 8: Phase evolution, (f>k(t), at wavenumber k = (25,0). 



each wavenumber k there are two modes: 

a k = ci + c 2 e~ i<JJ *\ 

where u* > u>k and the complex amplitudes are c\ and c 2 varying in time slower than the linear 
oscillations. Most of the time \c\\ > | C2 1 and, therefore, the phase <pk oscillates periodically about the 
some mean value (equal to the phase of ci). However, sometimes c 2 becomes greater than c\ and then 
the path of a k will encircle zero in the complex plane, so that (fik starts gaining 2n for each rotation. 
These are the phase runs, and in order to trigger these runs the complex-plane path of a k must encircle 
zero, which explains the observed small values of this quantity during the phase runs. This effect is 
related to the general phenomenon of phase singularities of complex fields at points of zero amplitude 
(e.g. [34]). However, our phase runs continue all the time until c 2 becomes greater than C\ and during 
this time interval (excluding its ends) the amplitude is non-zero and the phase is perfectly regular. 

Note that such observed behavior of the phase is very dependent on the wave variables: it occurs 
in natural variables such as the surface height and velocity but is absent for the waveaction variable 
after the canonical transformation leading to the Zakharov equation. Thus, the phase runs can say 
nothing about the dynamics or statistics of the waveaction amplitude used in WT, but they mark short 
time intervals when nonlinearity fails to be weak for a particular /c-mode, and this will be used as a 
diagnostic tool in our simulations. 

In figure 11 we show comparison of the mean rate of the phase change during the upward phase 
runs and the second peak's frequency at different wavenumbers. A significant coincidence of these two 
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Figure 9: Detailed graph of A\(t) at wavenumber k = (25, 0) corresponding to a time interval character- 
ized by low amplitudes and, therefore, phase runs. The linear wave period for this mode is 2ix /u k 1.256 
and the characteristic nonlinear time is of the same order of magnitude. 
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Figure 11: Mean rate of the phase change during the upward phase runs (thick curve) and the second 
peak's frequency (thin curve) as functions of wavenumber. 



curves supports the proposed above two-mode explanation. 

A word of caution is due about the simple two-mode explanation of the phase runs. Indeed, according 
to this picture the phase should always run to higher values whereas in figure 8 we see both upward and 
downward runs. A possible explanation of this is that phase runs may be triggered not only by sharp 
peaks but also by broadband distributions with frequencies less than the linear one. Such broadband 
distributions could be made of modes which are strong (and therefore can produce phase runs) but 
whose duration in time is short and sporadic and at different frequencies (hence a broad spectrum). 
This picture is supported by the fact that the amount of total wave energy in the frequency range below 
the linear frequency is similar (and for low wavenumbers even greater) than the amount of energy in 
the range above the linear frequencies, see figure 13. Evidence of relation of the sub-linear modes and 
the downward phase runs is shown in figure 12 where the mean rate of the phase change during the 
downward phase runs and the frequency of the highest sub-linear peak. Again, we can see agreement 
of these two curves although with a greater level of fluctuations due to the fact the sub-linear modes 
are very spread over different frequencies. 

7.5 Nonlinear ly active modes and cascade "avalanches" 

Component with the shifted frequency u* is clearly a nonlinear effect (there is no frequency shift in 
linear dynamics). Thus, the relative strength of c 2 and c\ can be used as a measure of nonlinearity. 
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Figure 12: Mean rate of the phase change during the downward phase runs (thick curve) and the 
frequency of the highest sub-linear peak (thin curve) as functions of wavenumber. 



Particularly, the phase runs mark the events when nonlinearity becomes strong. Figure 10 shows 
locations of the phase runs in the 2D wavenumber space which happened at t — 500. Note that at 
that time ZF steady spectrum has already formed. In the energy cascade range, we see that the phase 
run density is increasing toward high fc's, which is in agreement with the WT prediction that the 
nonlinearity grows as one cascades down-scale [9,10]. Curiously, we also observe high density of the 
phase runs within th circle k < 6, which is, perhaps, manifestation of a waveaction accumulation via 
an inverse cascade process. However, this range is too small for any meaningful conclusions to be made 
about the inverse cascade properties. 

The energy cascade from the forcing region toward the high wavenumber region proceeds in a non- 
uniform in time fashion somewhat resembling sporadic sandpile avalanches. This arises due to the 
/c-grid discreteness effects which tend to block the resonant wave interaction when the wave intensities 
are small. This situation resembles "frozen turbulence" of [14]. Thus, the wave energy does not cascade 
to high wavenumbers and it tends to accumulate near the forcing scales until the wave intensity is strong 
enough to restore the resonant interaction via the nonlinear resonance broadening. At this moment the 
energy cascade toward high wavenumbers sets in, and this leads to depletion of energy at the forcing 
scale, - "sandpile tips over" . In turn, depletion of energy at the forcing scale leads to blocking of the 
energy cascade, and the process continues in a repetitive manner. As a result system oscillates between 
the state of "frozen turbulence" and the state of "avalanche cascade". This behavior is illustrated in 
figure 14 which shows percentage of modes experiencing phase runs in two different wavenumber ranges 
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Figure 13: Total energy of modes with frequencies below the linear frequency (solid line) and the modes 
with frequencies above the linear frequency (dashed line). 
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Figure 14: Percentage of modes experiencing phase runs in the range 13 < \k\ < 29 (lower curve) and 
in the range 30 < \k\ < 45 (upper curve). 



13 < k < 29 and 30 < k < 45. One can see that the shapes of these two curves bear a great degree 
of similarity up to a certain time delay and a vertical shift in the second curve with respect to the 
first one. The vertical shift reflects the fact that the energy cascade gets stronger as it proceeds to 
large wavenumbers. The time delay, on the other hand indicates the direction and the character of 
the sporadic energy cascade. It shows that a higher (lower) nonlinear activity at low fc's after a finite 
delay causes a higher (lower) activity at higher fc's, which could be compared with propagation of an 
avalanche (quenching) down a sandpile. 



7.6 Correlations of phases, phase factors and amplitudes. 

WT closure relies on the RPA properties of the wave fields, i.e. that the amplitudes A k and the phase 
factors tpk are statistically independent variables. On the other hand, WT calculation for the phases 4>k 
shows that these quantities get correlated. In order to check these properties and predictions numeri- 
cally, let us introduce a function that measures the degree of statistical dependence (or independence) 
of some Fourier-space variables X(ki) and Y(k 2 ) : 

c , k k) (x(k 1 )r(k 2 ))-(x(k 1 ))(r(k 2 )) 

vWki)> - (X(k 1 ))V(^ 2 (k 2 )) - <nk 2 )> 2 ' 
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Figure 15: Two-point auto-correlations for the phases C^(ki, k) (dashed line) and for the phase factors 
C^(kx,k) (solid line) with one point fixed at k = (15,0). 



For example, we can examine to what degree amplitudes A and independent of the phase factors ip by 
looking at the function C y i i ^(k 1 ,k 2 ) for different values of k x and k 2 . Independence of the amplitudes 
at different wavenumbers can be examined by the auto-correlation function ^ ^(ki, k 2 ), and similar 
for the phase factors and the phases. We restrict ourselves with choosing ki = (15,0) and k 2 = (k,0) 
with k G (10,64). Figure 15 shows the values of correlators C<^(ki, k 2 ) and C^(ki,k 2 ) as functions 
of k. In agreement with WT predictions, auto-correlations of -?/Vs are very small whereas the ones of 
fc 's are significant (except, of course, for k = 15 where by definition these correlators are equal to one). 
Correlators (^^(ki, k 2 ) and C^^kx, k 2 ) are shown in figure 16. Again, we see a good agreement with 
the WT prediction: these correlations are very small (except, again, Ca,a(15, 15) = 1). 



8 Discussions 

In this paper, we used direct numerical simulations of the free water surface in order to examine 
the statistical properties of the water-wave field beyond the energy spectrum. Our first aim was to 
check recent predictions of the WT theory about the PDF and intermittency, about the character of 
correlations of the wave amplitudes and phases. We particularly focused on the question how the effects 
of discreteness and finite nonlinearity change statistics with respect to the WT closure developed for 
weak nonlinearities and for a continuous wavenumber space. 
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Figure 16: Two-point auto-correlation of amplitudes (^^(k^k) (thin curve) and two-point correlation 
between the amplitudes and phase factors Ca^O^i, k) (thick curve) with one point fixed at k = (15, 0). 



Firstly, following [1-3] we see formation of a quasi-steady spectrum consistent with the Zakharov- 
Filonenko spectrum predicted by WT. Secondly, we measured PDF for the wave amplitudes and observe 
an anomalously large, with respect to Gaussian fields, probability of strong waves. This result is 
in agreement with recent theoretical predictions of [4, 5] . Thirdly, we measure correlations for the 
amplitudes, phases and phase factors and we observe agreement with predictions of [4, 5] . Namely, the 
amplitude and the phase correlations behave as statistically independent variables, whereas the phases 
develop strong auto-correlations over the nonlinear time. Note that these properties are fundamental 
for the WT closure to work, so in a way we provide a numerical validation for the WT approach. 

We also find that at each k there are two sharp frequency peaks: a dominant one at the linear 
frequency and a weaker one with a frequency shift arising due to the k/2-mode. Somewhat related to 
this two-peak frequency structure is the observed time behavior of the phase. We observe calm periods 
during which the phase oscillates within 27r-wide margins intermittent with sudden phase runs during 
which it experiences a monotonic change significantly greater than 2tt. 

Finally, we observe that the energy cascade is "bursty" in time and is somewhat similar to sporadic 
sandpile avalanches. We give a plausible explanation of this behavior as an interplay of effects of 
discreteness and nonlinearity. Because in between of the avalanche discharges the resonances are absent 
then, at least qualitatively, one can refer to the KAM theory and say that the evolution should remain 
close to the corresponding integrable case, - the linear system in our case. This picture is supported 
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by a simple analysis of quasi-resonances given in this paper which indicates that there exists a single 
threshold value of turbulence intensity at the forcing scale separating the no-cascade and unlimited (in 
k) cascade regimes. 

A further numerical study of the avalanche effect is desirable, particularly using a different wavenum- 
ber grid and using a more direct method of measuring the turbulent flux and its correlations for different 
inertial interval points. 
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